

#'@export
#'@method print acifit

print.acifit <- function(x, include = c("data","rmse","parameters","method","citransition",
                                        "settings","gstar"), ...){

  cat("Result of fitaci.\n\n")

  if("data" %in% include){
    cat("Data and predictions:\n")
    print(x$df)
  }

  if("rmse" %in% include){
    cat("\nRoot mean squared error: ", x$RMSE, "\n")
  }

  if("parameters" %in% include){
    cat("\nEstimated parameters:\n")

    print(x$pars)
    if(x$Tcorrect)
      cat("Note: Vcmax, Jmax are at 25C, Rd is at measurement T.\n")
    else
      cat("Note: Vcmax, Jmax, Rd are at measurement T.\n")

    if(!is.na(x$gmeso)){
      cat("Note: Mesophyll conductance was input, Vcmax and Jmax are Cc-based rates.\n")
    }

    if(x$Rd_measured)
      cat("Note: measured Rd was provided, only Vcmax and Jmax were fit.\n")
  }

  if("method" %in% include){
    cat("\nCurve was fit using method: ", x$fitmethod, "\n")
  }


  if("citransition" %in% include){
    if(!is.na(x$citransition)){
      cat("\nCi transition was constrained to be: ", x$citransition, "\n")
      cat("Actual fitted Ci transition: ", x$Ci_transition,"\n")
    }
  }

  if("settings" %in% include){
    cat("\nParameter settings:\n")
    fm <- formals(x$Photosyn)
    pars <- c("Patm","alpha","theta","EaV","EdVC","delsC","EaJ","EdVJ","delsJ")
    fm <- unlist(fm[pars])
    cat(paste0(pars," = ", fm,"\n"))
  }

  if("gstar" %in% include){
    if(!x$gstarinput | !x$kminput){
      cat("\nEstimated from Tleaf (shown at mean Tleaf):\n")
      if(!x$gstarinput)cat("GammaStar = ",x$GammaStar,"\n")
      if(!x$kminput)cat("Km = ",x$Km,"\n")
    }

    if(x$gstarinput | x$kminput){
      cat("\nSet by user:\n")
      if(x$gstarinput)cat("GammaStar = ",x$GammaStar,"\n")
      if(x$kminput)cat("Km = ",x$Km,"\n")
    }
  }

}


#' @method summary acifit
#' @export
summary.acifit <- function(object,...){

  print.acifit(object, ...)

}



#' @method coef acifit
#' @export
coef.acifit <- function(object, ...){
  v <- unname(object$pars[,1])
  names(v) <- rownames(object$pars)
  return(v)
}






#' @method fitted acifit
#' @export
fitted.acifit <- function(object,...){

  object$df$Amodel

}



#' @method plot acifit
#' @export
#' @param x For plot.acifit, an object returned by \code{fitaci}
#' @param xlim Limits for the X axis, if left blank estimated from data
#' @param ylim Limits for the Y axis, if left blank estimated from data
#' @param whichA By default all photosynthetic rates are plotted (Aj=Jmax-limited
#' (blue), Ac=Vcmax-limited (red), Hyperbolic minimum (black)), TPU-limited rate
#' (Ap, if estimated in the fit). Or, specify one or two of them.
#' @param what The default is to plot both the data and the model fit, or specify 'data' or
#' 'model' to plot one of them, or 'none' for neither (only the plot region is set up)
#' @param add If TRUE, adds to the current plot
#' @param pch The plotting symbol for the data
#' @param addzeroline If TRUE, the default, adds a dashed line at y=0
#' @param addlegend If TRUE, adds a legend (by default does not add a legend if add=TRUE)
#' @param legendbty Box type for the legend, passed to argument bty in \code{\link{legend}}.
#' @param transitionpoint For plot.acifit, whether to plot a symbol at the transition point.
#' @param linecols Vector of three colours for the lines (limiting rate, Ac, Aj), if one value
#' provided it is used for all three.
#' @param lwd Line widths, can be a vector of length 2 (first element for Ac and Aj, second one
#' for the limiting rate).
#' @param lty Line type (only for Amin - the limiting rate).
#' @rdname fitaci
#' @importFrom graphics points
#' @importFrom graphics abline
#' @importFrom graphics legend
#' @importFrom graphics text
plot.acifit <- function(x, what=c("data","model","none"), xlim=NULL, ylim=NULL,
                        whichA=c("Ac","Aj","Amin","Ap"), add=FALSE, pch=19,
                        addzeroline=TRUE, addlegend=!add, legendbty='n',
                        transitionpoint=TRUE,
                        linecols=c("#E41A1C", "#377EB8", "#4DAF4A"),
                        lwd=c(2.5,2.5), lty=1,
                        ...){

  # Note that Ci on the X-axis is in molar units!
  if(is.null(ylim))ylim <- with(x$df, c(min(Ameas), 1.1*max(Ameas)))
  if(is.null(xlim))xlim <- with(x$df,c(0, max(Ci_original)))
  if(length(lwd)==1)lwd <- c(lwd,lwd)
  if(length(linecols)==1)linecols <- rep(linecols,3)

  # Vector of Ci values at which to evaluate fitted ACi curve.
  Ci <- with(x$df, seq(min(Ci_original), max(Ci_original), length=101))

  # Exact model used to fit the A-Ci curve was saved in the object.
  # (parameter settings etc. are preserved)
  pcor <- mean(x$df$Patm)/100

  pred <- x$Photosyn(Ci=Ci * pcor)
  pred$Ci_original <- pred$Ci / pcor

  # If TPU could not be estimated, Ap will be NA for Ci > 400.
  pred$Ap[is.na(pred$Ap)] <- 1000

  # Is there a TPU limitation?
  TPUlimit <- any(pred$Ap < pred$Aj)

  if(!add){
    with(x$df, plot(Ci_original, Ameas, type='n',
                    ylim=ylim,
                    xlim=xlim,
                    xlab=expression(italic(C)[i]~~(ppm)),
                    ylab=expression(italic(A)[n]~~(mu*mol~m^-2~s^-1)),
                    ...
    ))
  }
  if("data" %in% what)with(x$df, points(Ci_original, Ameas, pch=pch,...))

  if("model" %in% what){
    if("Aj" %in% whichA)with(pred, lines(Ci_original, Aj-Rd, col=linecols[2],lwd=lwd[1]))
    if("Ac" %in% whichA)with(pred, lines(Ci_original, Ac-Rd, col=linecols[3],lwd=lwd[1]))
    if("Ap" %in% whichA & TPUlimit){
      predp <- pred[pred$Ci_original > 400,]
      with(predp, lines(Ci_original, Ap-Rd, col="darkgrey", lty=5, lwd=lwd[1]))
    }
    if("Amin" %in% whichA)with(pred, lines(Ci_original, ALEAF, col=linecols[1], lwd=lwd[2], lty=lty))
  }

  if(transitionpoint && "model" %in% what){
    points(x$Ci_transition / pcor, x$Photosyn(Ci=x$Ci_transition)$ALEAF,
           pch=19, col="#984EA3", cex=1.1)

    points((x$Ci_transition / pcor) * 0.17, (x$Photosyn(Ci=x$Ci_transition)$ALEAF) * 1.1,
           pch=19, col="#984EA3", cex=0.8)

    text(
      (x$Ci_transition / pcor) * 0.18,
      (x$Photosyn(Ci = x$Ci_transition)$ALEAF) * 1.1,
      paste0("Ci transition = ", round(x$Ci_transition / pcor, 2)),
      cex = 0.6, pos = 4
    )
    }

  if(addzeroline){
    abline(h=0, lty=3)
  }

  if(addlegend & ! TPUlimit){
    legend("bottomright", c(expression(italic(A)[c]),
                            expression(italic(A)[j]),
                            "Limiting rate"), lty=1, lwd=c(lwd[1],lwd[1],lwd[2]),
          col= linecols[3:1], bty=legendbty, bg="white", cex = 0.7)
  }
  if(addlegend & TPUlimit){
    legend("bottomright", c(expression(italic(A)[c]),
                            expression(italic(A)[j]),
                            expression(italic(A)[p]),
                            "Limiting rate"
                            ), lty=1, lwd=c(rep(lwd[1],3), lwd[2]),
           col=c(linecols[3:2],"darkgrey",linecols[1]), bty=legendbty, bg="white",cex = 0.7)
  }

}
